 function [p_cal, T_cal,record_time]=THERMAID(frac_info,apert)
%THERMAID
%  data         ->  seismic data stored in txt format (optional, default 'inputFile')
%  sseis        ->  parameter for fracture generation 
%  attach_pre   ->  must be a string (optional, default 'attach_pre_timestep')
%  p_cal        ->  presure calculation at observation position
%  T_cal        ->  temperature calculation at observation positions

%% DEFAULTS AND GLOBAL PARAMETERS                                        %%
% global innerIter
showPlot = 1; 
attach_post = 'attach_post_timestep'; 

h=waitbar(0,'calculating');

%%        CREATE INPUT FILES
run('Matrix_Config')
run('Frac_Discrete')
run('Frac_Config')

%% define observation positions 
ntim=round(udata.timeSim/udata.dt)  ; 
npos= size(udata.obs, 1);
p_cal=zeros (ntim, npos-2);
T_cal=zeros (ntim, npos-2);
record_time =zeros(ntim,1);
%% ------------------------------------------------------------
x = linspace(0,udata.len(1),2*udata.Nf(1)+1); x = x(2:2:end);              % x grid vector 
y = linspace(0,udata.len(2),2*udata.Nf(2)+1); y = y(2:2:end);              % y grid vector

tNew  = udata.T0  .* ones(udata.Nf);                                       % initialize matrix   temperature solution vector
tNewf = udata.T0f .* ones(udata.Nf_f,1);                                   % initialize fracture temperature solution vector
tN    = [tNew(:); tNewf(:)];

p  = udata.p0  .* ones(udata.Nf);                                          % initialize matrix   pressure saturation solution vector
pf = udata.p0f .* ones(udata.Nf_f,1);                                      % initialize fracture pressure saturation solution vector
pN    = [p(:); pf(:)];

%%
CI = intersectionsGrid(udata,x,y,XY1);                                      % Find the intersection of each segment with the grid and compute the
[row,col] = intersectionsSegments(XY1,XY1);                                 % Connectivity Index based on the segments lengths and the average distance
                                                                            % to the fracture
%%
if ~(isempty(udata.ibcp))
    ind_ij = sub2ind(udata.Nf,udata.ibcp(:,1),udata.ibcp(:,2));            % In case of pressure controlled well BC find the fracture segments that
    for i = 1:length(ind_ij)                                               % lie in the given matrix cell where the BC is applied to
        ind = find(CI(:,1) == ind_ij(i),1);
        if ~(isempty(ind))
            udata.ibcp(i,3) = CI(ind,2);
        end
    end
end

%% ---------------- SIMULATION  --------------------------------- %%
time = udata.dt;                                                           % Set the time for the first timestep                                             
dt0 = udata.dt;
i = 0;                                                                     % Count all timesteps
while (time <= udata.timeSim)                                              % Time loop
    
%     if showPlot 
%         fprintf(['Simulation time = ' time2str(time)  '\n']);
     waitbar(time/udata.timeSim, h)
%     end

    tOld      = tNew;                                                      % Update saturation of previous iteration (time loop)
    tOldf     = tNewf;
    pOld      = p;                                                         % Update saturation of previous iteration (time loop)
    pOldf     = pf;
    epsP       = inf;
    epsT       = inf;
    i=i+1;
    innerIter = 1;
    while ((epsP >= udata.tol || epsT >= udata.tol) && innerIter <= udata.maxit) % Inner loop due to saturation dependence of gravity and viscosity
        tIt = tN;                                                          % Update saturation of previous iteration (inner loop)
        pIt = pN;
        
        try
            [Tx,Ty,Tf,Tfm,Tff,DfmT,g,gf,density_l,density_lf]  = initialize(udata,tNew,tNewf,p,pf,CI,row,col);
        catch
            warning('Matrix match error occured again!')
            return
        end                                                               % Initialize (update) the transmissivity based on new temperature
                                                                          
        [pN]= pressureSystem(udata,pOld,pOldf,Tx,Ty,Tf,g,gf,Q,Tfm,Tff);    % Construct pressure matrix and rhs vectors

        p = reshape(pN(1:prod(udata.Nf)),udata.Nf(1),udata.Nf(2));         % Assign solution to matrix solution array 
        pf = pN(prod(udata.Nf)+1:length(pN));                              % Assign solution to fracture solution vector

        [vx,vy,vf,Vfm,Vmf,Vff] = calcVelocity(udata,p,pf,Tx,Ty,Tf,Tfm,Tff,g,gf);  % Calculate Darcy velocities     
        
        if(udata.flagHeatTransport)
            [tN]  = transport_heat_System(udata,tOld,tOldf,vx,vy,vf,Vfm,Vmf,Vff,DfmT,Q,QT,density_l,density_lf); % Solve mass transport equation 

            tNew = tN(1:prod(udata.Nf));                     
            tNew= reshape(tNew,udata.Nf(1),udata.Nf(2));                   % Assign transport solution to matrix solution array 
            tNewf = tN(prod(udata.Nf)+1:length(tN));                       % Assign transport solution to fracture solution array 
        end

        epsP = norm((abs(pN(:) - pIt(:))),inf); 
        epsT = norm((abs(tN(:) - tIt(:))),inf); 

%         if (innerIter < 2) 
%             epsP0 = epsP; 
%         end
        
%         if showPlot
%             %fprintf('\t Residual at %d. loop: %d and %d\n', innerIter,epsP, epsT);   
%             fprintf('\t Residual at %d. \n', innerIter); 
%         end
        innerIter = innerIter+1;

        if (innerIter == udata.maxit)
            error('outer finescale loop did not converge'),
        end
    end
    
    p_cal(i,:)= p(udata.Nf(1)*(obs_xy(3:npos,2)-1)+obs_xy(3:npos,1))';
    T_cal(i,:)= tNew(udata.Nf(1)*(obs_xy(3:npos,2)-1)+obs_xy(3:npos,1))';
    record_time (i)= time+udata.dt;
    
 %    run (attach_post)

    %%     Time Step Control                                             %%

    if (innerIter > 25)
        udata.dt = 0.6*udata.dt;
%     elseif (innerIter < 4) % REMOVED THIS FOR RENO SIMULATIONS
%         udata.dt = min(1.2*udata.dt,dt0);
    end
        
    if (time+udata.dt > udata.timeSim) 
        udata.dt = udata.timeSim - time;
        if udata.dt == 0 
            udata.dt = nan;
        end
    end

    time=time+udata.dt;
end
delete(h)

% Write THERMAID.m  variables into the workspace after the simulation has ended. 

% ListOfVariables = who;
% for k = 1:length(ListOfVariables)
%    assignin('base',ListOfVariables{k},eval(ListOfVariables{k}))
% end
% 
% try
%     close(vidObj);
% catch
%      warning('No video object to close')

end
